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Abstract.  An  image  registration  technique  is  presented  for  the  registration 
of  medical  images  using  a  hybrid  combination  of  coarse-scale  landmark  and 
B-splines  deformable  registration  techniques.  The  technique  is  particularly  ef¬ 
fective  for  registration  problems  in  which  the  images  to  be  registered  contain 
large  localized  deformations.  A  brief  overview  of  landmark  and  deformable  reg¬ 
istration  techniques  is  presented.  The  hierarchical  multiscale  image  decompo¬ 
sition  of  E.  Tadmor,  S.  Nezzar,  and  L.  Vese,  A  multiscale  image  representation 
using  hierarchical  (BV,  L2)  decompositions ,  Multiscale  Modeling  and  Simula¬ 
tions,  vol.  2,  no.  4,  pp.  554-579,  2004,  is  reviewed,  and  an  image  registration 
algorithm  is  developed  based  on  combining  the  multiscale  decomposition  with 
landmark  and  deformable  techniques.  Successful  registration  of  medical  im¬ 
ages  is  achieved  by  first  obtaining  a  hierarchical  multiscale  decomposition  of 
the  images  and  then  using  landmark-based  registration  to  register  the  resulting 
coarse  scales.  Corresponding  bony  structure  landmarks  are  easily  identified  in 
the  coarse  scales,  which  contain  only  the  large  shapes  and  main  features  of  the 
image.  This  registration  is  then  fine  tuned  by  using  the  resulting  transforma¬ 
tion  as  the  starting  point  to  deformably  register  the  original  images  with  each 
other  using  an  iterated  multiscale  B-splines  deformable  registration  technique. 
The  accuracy  and  efficiency  of  the  hybrid  technique  is  demonstrated  with  sev¬ 
eral  image  registration  case  studies  in  two  and  three  dimensions.  Additionally, 
the  hybrid  technique  is  shown  to  be  very  robust  with  respect  to  the  location 
of  landmarks  and  presence  of  noise. 


1.  Introduction.  Image  registration  is  the  process  of  determining  the  optimal  spa¬ 
tial  transformation  that  brings  two  images  into  alignment  with  each  other.  Image 
registration  is  necessary,  for  example,  when  images  of  the  same  object  are  taken  at 
different  times,  from  different  imaging  devices,  or  from  different  perspectives.  The 
two  images  to  be  registered,  called  the  fixed  and  moving  images,  are  the  input  to 
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the  registration  algorithm;  the  output  is  the  optimal  transformation  that  maps  the 
moving  image  to  the  fixed  image.  Ideally,  the  transformed  moving  image  should  be 
identical  to  the  fixed  image  after  registration.  Applications  of  image  registration 
include  image-guided  radiation  therapy  (IGRT),  image-guided  surgery,  functional 
MRI  analysis,  and  tumor  detection,  as  well  as  many  non-medical  applications,  such 
as  computer  vision,  pattern  recognition,  and  remotely  sensed  data  processing  (see 
[5]  and  the  references  therein). 

Image  registration  models  are  classified  into  two  main  categories  according  to  the 
transformation  type:  rigid  and  deformable.  Rigid  image  registration  models  assume 
that  the  transformation  that  maps  the  moving  image  to  the  fixed  image  consists  only 
of  translations  and  rotations,  while  deformable  models  allow  localized  stretching 
of  images.  Rigid  models  are  sufficient  in  certain  circumstances.  However,  many 
registration  problems,  particularly  in  medical  applications,  are  non-rigid  since  most 
of  the  organs  in  the  human  body  are  not  confined  to  rigid  motion.  For  example, 
respiratory  motion  causes  non-rigid,  or  deformable,  distortion  of  the  lungs  and 
other  organs.  As  another  example,  image-guided  neurosurgery  procedures  require 
deformable  registration  of  pre-  and  intra-operative  images  of  the  brain  [18],  [26]. 
For  additional  applications  of  deformable  registration,  see  [11],  [22],  [23],  [24],  and 
the  references  therein.  Most  current  research  in  image  registration  is  focused  on  the 
problem  of  deformable  registration,  and  that  is  our  focus  in  this  paper. 

We  present  a  novel  algorithm  for  the  deformable  registration  of  medical  images 
that  combines  landmark  and  deformable  registration  models  using  a  hierarchical 
multiscale  decomposition  of  the  images  to  be  registered.  Landmark-based  registra¬ 
tion  techniques  use  the  correspondence  of  a  set  of  features,  or  landmarks,  in  the 
images  to  determine  the  transformation  that  maps  the  moving  image  to  the  fixed 
image.  Although  landmark-based  techniques  are  computationally  easy  to  imple¬ 
ment,  the  identification  of  corresponding  features  in  the  images  to  be  registered  is  a 
difficult  and  time-consuming  task,  and  the  accuracy  of  such  techniques  is  dependent 
on  precise  correspondence  between  landmarks  [20].  Our  hybrid  technique,  however, 
uses  a  combination  of  coarse-scale  landmark-based  and  deformable  registration  tech¬ 
niques,  and  is  not  dependent  upon  exact  correspondence  of  landmarks.  The  hybrid 
technique  is  shown  to  be  particularly  robust  to  perturbations  in  the  location  of 
landmarks,  and  thus  is  a  significant  improvement  over  ordinary  landmark-based 
registration  techniques.  The  development  of  our  technique  is  motivated  by  the  idea 
that  correspondence  between  certain  structures  in  images,  such  as  bony  structures, 
can  be  easily  identified  visually,  while  the  correspondence  between  other  regions, 
such  as  soft  tissue,  is  less  easy  to  see.  Bony  structures  undergo  rigid  transfor¬ 
mations;  soft  tissue  and  muscular  structures  undergo  deformable  transformations. 
Thus,  mapping  of  the  two  regions  should  be  approached  using  different  techniques. 

This  paper  is  an  extension  of  [15]  and  [16],  in  which  we  presented  a  multiscale 
approach  to  the  problems  of  rigid  and  deformable  registration  of  images  that  con¬ 
tain  large  levels  of  noise.  This  paper  is  organized  in  sections.  In  Section  2,  we 
provide  a  brief  overview  of  the  image  registration  problem  and  discuss  landmark 
and  deformable  registration  techniques.  In  Section  3,  we  review  the  hierarchical 
multiscale  image  decomposition  of  [25],  and  in  Section  4  we  develop  a  hybrid  im¬ 
age  registration  algorithm  based  on  combining  the  decomposition  with  landmark 
and  deformable  registration  techniques.  The  accuracy  and  efficiency  of  our  hybrid 
technique  is  demonstrated  in  Section  5  with  several  image  registration  case  studies 
in  both  two  and  three  dimensions.  In  Section  6,  we  demonstrate  the  robustness  of 
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our  algorithm  with  respect  to  the  location  of  landmarks,  number  of  landmarks,  and 
presence  of  noise.  Concluding  remarks  are  given  in  Section  7. 

2.  The  registration  problem.  Image  registration  is  the  problem  of  determining 
the  optimal  spatial  transformation  between  two  given  images.  The  images  to  be 
registered  are  typically  (though  not  necessarily)  images  of  the  same  object  acquired 
at  different  times,  from  different  perspectives,  or  from  different  modalities.  Image 
registration  has  many  important  applications  in  medicine  (e.g.,  diagnosis,  analysis, 
and  treatment  planning/evaluation),  astrophysics  (e.g.,  alignment  of  images  from 
different  frequencies),  military  applications  (e.g.,  target  recognition),  computer  vi¬ 
sion,  remote  sensing,  and  many  others.  See  [8],  [13]  for  an  overview  of  applications 
of  image  registration.  In  this  paper,  we  are  primarily  concerned  with  registration 
of  medical  images. 

2.1.  The  mathematical  setting.  An  n-dimensional  gray-scale  image  /  is  a  map¬ 
ping  which  assigns  to  every  point  x  G  £2  C  Mn  a  gray  value  /(x),  called  the  intensity 
value  of  the  image  /  at  the  point  x.  Given  two  images  A(x)  and  B(x ),  referred  to  as 
the  fixed  and  moving  images,  the  goal  of  the  registration  problem  is  to  find  a  spatial 
transformation  <j)  such  that  the  fixed  image  A{x)  and  the  transformed  moving  image 
B((j)(x ))  are  similar  in  some  appropriate  mathematical  way. 

A  registration  algorithm  has  three  components:  the  transformation  model  which 
specifies  the  way  in  which  the  moving  image  can  be  transformed  to  correspond  to  the 
fixed  image,  the  similarity  metric  used  to  compare  the  images,  and  the  optimization 
process  that  varies  the  parameters  of  the  transformation  model  in  such  a  way  that 
the  resulting  transformation  is  optimal.  Given  a  distance  metric  D  on  L2(M2) 
(for  two-dimensional  images),  and  a  specified  space  of  transformation  models,  the 
general  registration  problem  can  be  stated  as  follows: 

<p  =  argmin  D(A ,  B( 

where  ip  is  in  the  specified  space  of  transformation  models.  Common  metrics  used  in 
image  registration  are  mean  squares,  normalized  correlation,  and  mutual  informa¬ 
tion  (for  registration  of  multi- modality  images).  The  space  of  transformation  models 
used  in  a  particular  registration  problem  depends  on  the  physical  and  anatomical 
characteristics  of  the  images  to  be  registered.  For  example,  rigid  and  affine  transfor¬ 
mations  are  used  when  the  moving  image  is  assumed  to  differ  from  the  fixed  image 
by  translation,  rotation,  dilation,  and  shear.  Polynomial  and  splines  translations 
are  used  when  the  transformation  between  the  images  is  assumed  to  consist  of  local¬ 
ized  stretching  and  non-rigid  deformation.  See  [5],  [13]  for  an  overview  of  the  image 
registration  problem.  In  this  paper,  we  will  use  a  combination  of  landmark-based 
and  deformable  registration  techniques. 

2.2.  Landmark-based  registration.  Landmark-based  registration  algorithms  use 
the  (manual  or  automatic)  identification  of  corresponding  anatomical  structures  (or 
other  features)  in  the  images  to  be  registered.  For  an  overview  of  landmark-based 
registration,  see  [20].  Given  two  images  A  and  B  to  be  registered  with  one  another, 
let  Fk(A)  and  Fk(B),  k  =  1,2,...  m  denote  m  corresponding  features  of  the  im¬ 
ages.  The  solution  <p  of  the  registration  problem  is  then  a  map  <p  :  M2  — >  M2  (in  the 
two-dimensional  case)  such  that 


Fk(A)  =  </)(Fk(B)),  k  =  l,...m. 


(1) 
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More  generally,  if  ||  •  ||  is  a  norm  on  the  features  space,  the  landmark-based  regis¬ 
tration  problem  can  be  stated  as  the  following  minimization  problem: 

4>  =  argmin  DLM  (-0),  (2) 

ijj:R2^R2 

m 

where  DLM  (ip)  :=  ^  \\Bk{A)  —  ^(F^B) ||2.  To  solve  the  minimization  given  by 
k= 1 

Eq.  (2),  the  transformations  *ip  are  typically  chosen  to  be  an  element  of  some  finite¬ 
dimensional  space,  such  as  polynomials,  splines,  or  wavelets,  and  the  minimization 
can  be  easily  solved  by  expanding  DLM  (<p>)  in  terms  of  basis  functions  of  the  space 
and  solving  the  resulting  least-squares  problem. 

Landmark-based  registration  techniques  are  computationally  simple  and  efficient, 
but  the  main  drawback  of  such  techniques  is  that  they  rely  on  precise  correspondence 
of  features  in  the  images  to  be  registered.  Although  there  has  been  recent  research 
on  automatic  identification  of  landmarks  (see  [2],  [7],  [24]),  in  practice  landmarks 
are  typically  identified  manually.  Precise  identification  of  corresponding  landmarks 
is  time  consuming  and  tedious,  even  for  a  medical  expert  [20].  In  addition,  there  are 
numerous  known  examples  of  cases  in  which  the  result  of  the  landmark  registration 
process  is  a  transformation  which  correctly  matches  the  user-supplied  landmarks 
but  is  not  physically  meaningful.  See  [14,  p.  44]  for  a  simple  example  of  such  a 
case. 

2.3.  Deformable  registration.  There  are  many  existing  approaches  to  deformable 
registration:  [1],  [4],  [12],  [19],  [23].  Splines-based  transformation  models  are  among 
the  most  common  and  important  transformation  models  used  in  non-rigid  registra¬ 
tion  problems  [6].  Splines-based  registration  algorithms  use  control  points  in  the 
fixed  and  moving  images  and  a  splines  function  to  define  transformations  away  from 
these  points.  The  two  main  spline  models  used  in  image  registration  problems  are 
thin-plate  splines  and  B-splines.  Thin-plate  splines  have  the  property  that  each 
control  point  has  a  global  influence  on  the  transformation.  That  is,  if  the  position 
of  one  control  point  is  perturbed,  then  all  other  points  in  the  image  are  perturbed 
as  well.  This  can  be  a  disadvantage  because  it  limits  the  ability  of  the  transfor¬ 
mation  model  to  model  localized  deformations.  In  addition,  the  computation  time 
required  for  a  thin-plate  spline-based  registration  algorithm  increases  significantly 
as  the  number  of  control  points  increases.  See  [3]  for  an  overview  of  thin-plate 
splines. 

In  contrast,  B-splines  are  defined  locally  in  the  neighborhood  of  each  control 
point.  Thus  perturbing  the  position  of  one  control  point  affects  the  transforma¬ 
tion  only  in  a  neighborhood  of  that  point.  As  a  result,  B-spline-based  registration 
techniques  are  more  computationally  efficient  than  thin-plate  splines,  especially  for 
a  large  number  of  control  points.  See  [9]  and  [10]  for  a  detailed  description  of  B- 
splines  transformation  models.  The  displacement  of  a  node  a^j  is  specified  by  a 
vector  x^j,  and  the  displacement  vectors  of  a  collection  of  nodes  characterize  the 
tissue  deformation.  The  displacement  at  a  particular  location  on  the  image  is  de¬ 
duced  by  fitting  a  polynomial,  expressed  using  the  basis  splines  (B-splines),  to  the 
grid  nodes.  The  number  of  control  points  used  determines  the  number  of  degrees 
of  freedom  of  the  transformation  model,  and  hence,  the  computational  complexity 
of  the  B-splines  algorithm. 

To  define  the  splines-based  deformation  model  in  two  dimensions,  we  let  Q  = 
{(x,y)  |  0  <  x  <  X,  0  <  y  <  Y}  denote  the  domain  of  the  image.  Let  a  denote 
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a  nx  x  ny  mesh  of  control  points  a^j  with  uniform  spacing  5.  Then  the  B-splines 
deformation  model  can  be  written  as  the  2-D  tensor  product  of  1-D  cubic  B-splines: 


2  2 


(j)(x,y)  =  EE  Bi  (x).Bm  (i/)c^_|_j 


(3) 


1=0  m= 0 


where  i  =  [x/nx\  —  1  ,  j  =  [y/ny\  —  1,  and  Bi  represents  the  Z-th  basis  of  the 
B-splines: 


B0(u)  =  hl-u)3  , 

6 

B\  (u)  =  ^  (3 u3  -  6 u2  +  4)  , 

6 

B2  (u)  =  —  ( — 3v?  T  3v?  T  3 u  T  1)  . 
6 


Changing  the  control  point  a^j  affects  the  transformation  only  in  a  local  neighbor¬ 
hood  of  aij.  The  control  points  a  act  as  parameters  of  the  B-splines  deformation 
model,  and  the  degree  of  non-rigid  deformation  which  can  be  modeled  depends  on 
the  resolution  of  the  mesh  of  control  points  a.  A  large  spacing  of  control  points 
allows  modeling  of  global  non-rigid  deformation,  while  a  small  spacing  of  control 
points  allows  modeling  of  local  non-rigid  deformations.  Additionally,  the  number  of 
control  points  determines  the  number  of  degrees  of  freedom  of  the  transformation 
model,  and  hence,  the  computational  complexity. 

The  splines  deformation  model  is  illustrated  schematically  in  Figure  1.  The 
deformation  is  computed  explicitly  at  control  points  by  computing  the  vector  dis¬ 
placement  between  the  fixed  and  moving  images,  and  the  splines  function  is  used 
to  interpolate  the  deformation  away  from  the  control  points.  The  implementation 
of  the  B-splines  deformation  registration  algorithm  works  in  the  following  way:  at 
each  iteration,  the  distance  D  between  the  fixed  and  moving  images  is  computed; 
the  parameters  of  the  B-splines  deformation  are  computed  based  on  the  direction 
of  steepest  gradient  descent  of  the  metric  D,  and  the  resulting  transformation  </>  is 
applied  to  the  moving  image;  the  distance  between  the  fixed  and  deformed  moving 
images  D(A ,  B((j)))  is  then  recomputed,  and  this  process  continues  until  the  distance 
D  between  the  images  is  optimized. 

3.  The  multiscale  decomposition.  The  multiscale  registration  techniques  to  be 
discussed  in  this  paper  are  based  on  the  multiscale  image  representation  using  the 
hierarchical  (BV,  L2)  decompositions  of  [25].  This  multiscale  decomposition  will 
provide  a  hierarchical  expansion  of  an  image  that  separates  the  essential  features 
of  the  image  (such  as  large  shapes  and  edges)  from  the  fine  scales  of  the  image 
(such  as  details  and  noise).  The  decomposition  is  hierarchical  in  the  sense  that 
it  will  produce  a  series  of  expansions  of  the  image  that  resolve  increasingly  finer 
scales,  and  hence  include  increasing  levels  of  detail.  The  mathematical  spaces  L2, 
the  space  of  square-integrable  functions,  and  BV ,  the  space  of  functions  of  bounded 
variation,  will  be  used  in  the  decomposition: 
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Figure  1.  Schematic  visualization  of  the  splines  deformation  model. 


BV={f  |  \\f\\BV- 

sup|h|_1||/(-  +  /i)  -  /(•) ||L1  <  oo}. 

0 

Generally,  images  can  be  thought  of  as  being  elements  of  the  space  L2(M2),  while 
the  main  features  of  an  image  (such  as  edges)  are  in  the  subspace  F?V(M2).  The 
multiscale  image  decomposition  of  [25]  interpolates  between  these  spaces,  providing 
a  decomposition  in  which  the  coarsest  scales  are  elements  of  BV  and  the  finest  scales 
are  elements  of  L2.  More  precisely,  the  decomposition  is  given  by  the  following. 
Define  the  J-functional  J(/,  A)  as  follows: 

A/A)  :=  _inf  AMiL  +  IMIw,  (4) 

where  A  >  0  is  a  scaling  parameter  that  separates  the  L2  and  BV  terms.  This 
functional  J(/,  A)  was  introduced  in  the  context  of  image  processing  by  Rudin, 
Osher,  and  Fatemi  [21].  Let  [u\,v\\  denote  the  minimizer  of  J(/,  A).  The  BV 
component,  u\,  captures  the  coarse  features  of  the  image  /,  while  the  L2  component, 
v\,  captures  the  finer  features  of  /  such  as  noise.  This  model  is  effective  in  denoising 
images  while  preserving  edges,  though  it  requires  prior  knowledge  on  the  noise 
scaling  A. 

Tadmor,  et  al.  proposed  in  [25]  an  alternative  point  of  view  in  which  the  mini¬ 
mization  of  J(/,  A)  is  interpreted  as  a  decomposition  /  =  u\  +  v\,  where  u\  extracts 
the  edges  of  /  and  v\  extracts  the  textures  of  /.  This  interpretation  depends  on  the 
scale  A,  since  texture  at  scale  A  consists  of  edges  when  viewed  under  a  refined  scale. 
We  refer  to  v\  =  f  —  u\  as  the  residual  of  the  decomposition.  Upon  decomposing 
f  =  u\  +  v a,  we  proceed  to  decompose  v\  as  follows: 


v\  =  u2\  +u2a, 


where 
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[m2a,v2a]  =  arginf  J(v\,  2A). 

U-\-V=V\ 

Thus  we  obtain  a  two-scale  representation  of  /  given  by  /  =  u\  +  R2A5  where  now 
v2\  —  f  —  (u\  +  ^2a)  is  the  residual.  Repeating  this  process  results  in  the  following 
hierarchical  multiscale  decomposition  of  /.  Starting  with  an  initial  scale  A  =  Ao, 
we  obtain  an  initial  decomposition  of  the  image  /: 

f  =  u0  +  v0,  [uo,  Vo]  =  arginf  J(/,  A0). 

U+V  =  f 

We  then  refine  this  decomposition  to  obtain 


Vj  =  uj+1  +  vj+i,  [Uj+i,vj+i\  =  arginf  J(vj,\02J+1),  j  =  0,1,... 

U-\-V=Vj 

After  k  steps  of  this  process,  we  have: 

/  =  Uq  +  Vq  =  Uq  +  U\  +  V\  =  Uq  +  U\  +  U2  +  V2  =  •  •  •  =  Uq  +  U\  +  .  .  .  +  Ufc  +  Vj^ ,  (5) 

which  is  a  multiscale  image  decomposition  /  ~  uq  +  u\  + . . .  with  a  residual  u/e. 
As  k  increases,  the  Uk  components  resolve  edges  with  increasing  scales  A&  =  Xo2k . 

4.  Hybrid  Landmark  and  deformable  registration  using  the  multiscale 
decomposition.  In  our  previous  work  [15],  [16],  we  used  the  multiscale  decom¬ 
position  of  [25]  reviewed  in  Section  3  to  develop  multiscale  rigid  and  deformable 
registration  algorithms  for  registration  of  images  that  contain  high  levels  of  noise. 
We  demonstrated  that  our  multiscale  registration  techniques  enable  successful  reg¬ 
istration  of  images  that  contain  noise  levels  significantly  greater  than  the  levels  at 
which  ordinary  rigid  and  deformable  registration  techniques  fail.  In  this  paper, 
we  extend  our  previous  work  on  multiscale  registration  to  a  new  registration  tech¬ 
nique  based  on  combining  landmark  and  deformable  registration  methods  with  the 
multiscale  decomposition.  We  shall  refer  to  this  technique  as  a  hybrid  multiscale 
landmark  and  deformable  registration  algorithm.  The  main  idea  behind  our  new 
technique  is  that  by  incorporating  landmark  registration,  we  can  use  known  infor¬ 
mation  about  correspondence  of  bony  structures  or  other  anatomical  structures  to 
improve  the  accuracy  and  efficiency  of  the  registration  procedure. 

For  the  general  setup,  consider  two  images  A  (the  fixed  image)  and  B  (the  moving 
image).  We  first  apply  the  multiscale  hierarchical  decomposition  to  both  images. 
Let  m  denote  the  number  of  hierarchical  steps  used  for  the  multiscale  decompo¬ 
sitions.  For  the  applications  presented  in  this  paper,  we  use  m  =  8  hierarchical 
scales.  Upon  decomposing  both  of  the  images  to  be  registered  we  use  the  coarse 
scales  of  each  image  to  identify  several  corresponding  anatomical  landmarks  in  each 
image.  Since  the  coarse  scales  contain  only  the  main  shapes  and  general  features 
of  each  image,  the  identification  of  landmarks  in  the  coarse  scales  is  a  relatively 
easy  task  and  can  potentially  be  implemented  automatically  (for  example,  using 
the  methods  of  [24]).  The  determination  of  which  coarse  scale  to  use  (e.g.,  the 
first  or  second  coarse  scale)  is  image  dependent;  upon  decomposing  each  image,  we 
choose  the  coarse  scale  that  contains  enough  general  shapes  to  identify  correspond¬ 
ing  landmarks  but  contains  few  details  and  noise.  For  the  medical  image  registration 
applications  presented  in  this  paper  it  is  sufficient  to  use  the  first  coarse  scale  of 
each  image. 
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To  identify  landmarks  we  must  choose  anatomical  features  of  the  two  images  that 
we  know  must  correspond  to  one  another.  For  example,  bony  structures  (such  as  the 
spine  or  ribs  in  images  of  the  torso)  in  one  image  must  correspond  to  bony  structures 
in  another  image.  We  then  use  the  identified  landmarks  to  perform  a  landmark- 
based  registration  of  the  coarse  scales  of  each  image.  The  output  of  the  landmark 
registration  algorithm  is  the  transformation  ^landmark  that  optimally  brings  the  first 
coarse  scales  C\(A)  and  C\(B)  into  spatial  alignment  with  one  another,  and  is  thus 
a  reasonable  approximation  for  the  optimal  transformation  between  the  original 
images  A  and  B.  Therefore,  we  use  landmark  as  the  starting  point  to  deformably 
register  the  next  coarse  scales  02(A)  and  02(B)  with  one  another.  We  then  use 
the  resulting  transformation  02  as  the  starting  point  to  register  the  next  scales 
03(A)  and  03(B) ,  and  iterate  this  procedure  until  the  final  scales  are  reached.  To 
summarize,  the  implementation  of  the  hybrid  multiscale  landmark  and  deformable 
registration  algorithm  is  as  follows: 


1.  Decompose  the  fixed  and  moving  images  to  be  registered  using  the  hierarchical 
multiscale  image  decomposition  of  [25]  reviewed  in  Section  3. 


Coarsest  Scale  of 
Image  A 


Finest  Scale  of 
Image  A 


Coarsest  Scale  of 
Image  B 


Finest  Scale  of 
Image  B 


2.  Register  the  coarse  scales  C\(A)  and  C\(B)  with  each  other  using  a  standard 
landmark-based  registration  algorithm.  This  step  allows  the  practitioner  to 
incorporate  known  anatomical  information  about  the  images  to  be  registered 
(such  as  correspondence  of  bony  structures)  into  the  registration  process.  Let 
^landmark  denote  the  resulting  transformation. 


Ci{A)  C\(B) 


•  1 

•3 4 

Landmark 

Registration 

•  1 

•  3 

•  2 

•  4 

0 landmark 

•  2 

•  4 

3.  Use  ^landmark  as  the  starting  point  to  deformably  register  the  next  scales  02(A) 
and  02(B)  with  each  other.  Let  0 2  denote  the  resulting  transformation. 

4.  Use  02  as  a  starting  point  to  deformably  register  03(A)  and  03(B)  with  each 
other.  Let  03  denote  the  transformation  obtained  upon  registering  03(A) 
with  03(B).  Iterate  this  method,  at  each  stage  using  the  transformation 
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computed  by  the  previous  scale  registration  algorithm  as  the  starting  point 
for  the  current  registration. 


0 landmark 

Ci(A) 

Deformable  ^ 

Registration 

Ci(B) 

01 

01 

Deformable 

c2(a) 

Registration 

Ci(B) 

02 

0m— 1 

Deformable 
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Registration 

Cm{B ) 
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5.  Results  and  discussion.  In  this  section,  we  demonstrate  the  accuracy,  effi¬ 
ciency,  and  robustness  of  the  image  registration  algorithm  presented  in  Section  4 
with  image  registration  experiments  using  both  clinical  and  synthetic  deformations. 
The  images  used  in  this  section  were  acquired  with  a  GE  Discover-ST  Scanner  (GE 
Medical  Systems,  Miluakee,  WI)  at  the  Stanford  University  Medical  Center.  We 
obtained  eight  sets  ( phase  bins )  of  images,  each  consisting  of  80  two-dimensional 
computed  tomography  (CT)  images  ( slices )  of  the  lungs.  Each  phase  corresponds 
to  a  different  breathing  phase  of  the  respiratory  cycle.  Each  2D  image  slice  con¬ 
tains  512  x  512  voxels,  and  the  slice  thickness  for  each  phase  is  2.5-mm,  and  the 
eight  breathing  phases  recorded  contain  approximately  400  MB  of  data  in  DICOM 
image  format.  Figure  2  illustrates  the  CT  image  slices  corresponding  to  the  first 
and  eighth  breathing  phases  of  the  respiratory  cycle. 

5.1.  2D  registration  of  respiratory  phases.  In  this  section,  we  demonstrate  the 
results  obtained  upon  registering  two  corresponding  sample  slices  (slice  16)  from  the 
first  and  eighth  phases  of  the  CT  data  set  using  both  ordinary  deformable  registra¬ 
tion  and  the  multiscale  hybrid  registration  algorithm.  The  slices  are  illustrated  in 
Figure  3.  For  ease  of  notation,  we  let  P1S16  denote  the  fixed  image  (on  the  left) 
and  P8S16  denote  the  moving  image  (on  the  right).  To  implement  the  multiscale 
hybrid  algorithm,  we  first  decompose  each  of  the  images  into  m  =  8  hierarchical 
scales,  and  identify  four  corresponding  pairs  of  landmarks  in  the  first  coarse  scales 
of  each  image,  as  shown  in  Figure  4. 

Figure  5  illustrates  the  difference  between  slices  before  and  after  registration  that 
we  performed  using  the  multiscale  hybrid  registration  technique.  For  comparison, 
we  also  illustrate  the  difference  between  the  slices  after  ordinary  deformable  regis¬ 
tration.  The  visual  results  presented  in  Figure  5  clearly  demonstrate  the  accuracy 
of  the  multiscale  hybrid  registration  technique.  The  hybrid  algorithm  successfully 
recovers  the  deformation  between  the  two  images,  and  is  more  accurate  than  ordi¬ 
nary  deformable  registration.  Moreover,  the  hybrid  algorithm  is  computationally 
more  efficient  than  ordinary  deformable  registration  techniques.  Working  on  a  Dell 
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Respiratory  Phase  1 


Respiratory  Phase  8 


Figure  2.  The  CT  image  slices  (80  per  phase)  corresponding  to 
the  first  and  eighth  breathing  phases  of  the  respiratory  cycle. 


Dimension  8400  Intel  Pentium  4  CPU  (3.40  GHz,  2.00  GB  of  RAM),  ordinary  de¬ 
formable  registration  of  the  images  shown  in  Figure  3  requires  approximately  115 
seconds.  The  total  time  required  for  the  hybrid  registration  algorithm  (including 
the  multiscale  decomposition  of  both  images,  landmark  registration  of  the  coarse 
scales,  and  the  final  deformable  registration)  is  approximately  72  seconds. 

To  further  quantitatively  evaluate  the  accuracy  of  the  hybrid  registration  tech¬ 
nique,  we  use  the  multiscale  hybrid  algorithm  to  register  all  80  slices  of  the  first 
breathing  phase  with  the  corresponding  slices  of  the  eighth  breathing  phase,  and  for 
each  registration  we  compute  the  correlation  coefficient  p  between  the  slices  before 
and  after  registration.  The  correlation  coefficient  p(A,B)  between  two  images  A 
and  B  is  given  by: 
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Figure  3.  Two  corresponding  sample  slices  from  two  breathing 
phases  of  the  same  patient.  The  image  on  the  left  (denoted  P1S16) 
is  from  the  first  breathing  phase,  and  the  image  on  the  right  (de¬ 
noted  P8S16)  is  the  corresponding  slice  from  the  eighth  breathing 
phase. 


Figure  4.  The  coarse-scale  landmarks  used  to  implement  the 
landmark-based  registration  portion  of  the  hybrid  algorithm. 


EE (Amn  -  A)(Bmn  -  B) 

P(A,B)  =  m  w 

WE  EOmn  -  A)2(Bmn  -  B)2 

y  ran 

where  A  and  B  are  mxn  two-dimensional  images  and  A  and  B  represent  the  mean 
value  of  the  elements  of  A  and  £?,  respectively.  A  correlation  coefficient  of  zero 
indicates  a  low  degree  of  matching  between  the  images,  and  a  correlation  coeffi¬ 
cient  of  1  indicates  exact  similarity  between  the  images.  Correlation  coefficients  are 
commonly  used  representations  of  similarity  between  images  for  the  evaluation  of 
deformable  registration  techniques  [17].  In  Figure  6,  we  plot  the  correlation  coef¬ 
ficients  between  the  slices  before  registration,  after  hybrid  multiscale  registration, 
after  multiscale  deformable  registration  (iteratively  registering  the  scales  of  the  im¬ 
ages  with  each  other  without  using  the  initial  landmark-based  registration)  [16], 
after  ordinary  deformable  registration,  and  after  landmark  registration. 

The  correlation  coefficients  plot  in  Figure  6  quantitatively  confirms  the  accuracy 
of  the  multiscale  hybrid  registration  technique.  For  all  80  CT  lung  slices  consid¬ 
ered,  the  hybrid  technique  significantly  improves  the  correlation  coefficient  between 
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Difference  Difference  After  Difference  After 

Before  Ordinary  Deformable  Hybrid  Multiscale 

Registration  Registration  Registration 


Figure  5.  The  difference  between  the  CT  respiratory  phase  slices 
P1S16  and  P8S16  before  registration,  after  ordinary  deformable 
registration,  and  after  hybrid  multiscale  registration. 


corresponding  images,  indicating  that  the  algorithm  successfully  recovers  the  defor¬ 
mation  between  the  breathing  phases.  Moreover,  the  hybrid  multiscale  technique 
is  more  accurate  (based  on  comparison  of  correlation  coefficients)  than  all  other 
registration  methods  considered  here. 

5.2.  Large  deformation  registration.  In  this  section,  we  demonstrate  the  ac¬ 
curacy  of  the  hybrid  technique  for  registration  images  that  contain  large  localized 
deformations.  We  apply  a  large  splines  deformation  to  the  fixed  image  P1S16  (slice 
16,  first  breathing  phase).  The  corresponding  deformation  field  and  the  fixed  and 
deformed  images  are  illustrated  in  Figure  7. 

As  in  Section  5.1,  we  first  decompose  the  images  to  be  registered  into  m  =  8 
hierarchical  scales.  Using  the  first  coarse  scale  of  each  image,  we  identify  four  corre¬ 
sponding  pairs  of  bony  structure  landmarks,  and  use  landmark-based  registration  to 
register  the  coarse  scales  of  the  images  with  one  another.  We  then  apply  the  result¬ 
ing  transformation  to  the  original  deformed  image,  and  use  a  B-splines  deformable 
registration  algorithm  to  complete  the  registration  process.  Figure  8  compares  the 
difference  between  the  images  before  registration,  after  ordinary  deformable  regis¬ 
tration  (for  comparison),  and  after  hybrid  multiscale  registration. 

The  results  presented  in  Figure  8  demonstrate  that  the  hybrid  registration  tech¬ 
nique  successfully  registers  images  that  contain  large  localized  deformations.  For 
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Figure  6.  The  correlation  coefficients  between  all  80  image  slices 
of  the  first  and  eighth  breathing  phase  before  registration,  after 
hybrid  multiscale  registration,  after  multiscale  deformable  regis¬ 
tration,  after  ordinary  deformable  registration,  and  after  landmark 
registration. 


Figure  7.  The  deformation  field  (left)  corresponding  to  the  de¬ 
formation  transformation  between  the  fixed  (center)  and  deformed 
(right)  images  to  be  registered  with  one  another. 
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Figure  8.  The  difference  between  the  fixed  and  deformed  images 
before  registration,  after  ordinary  deformable  registration,  and  af¬ 
ter  hybrid  multiscale  registration. 


such  images,  the  multiscale  hybrid  registration  technique  is  significantly  more  ac¬ 
curate  than  ordinary  deformable  registration.  The  improvement  of  the  multiscale 
hybrid  technique  over  ordinary  deformable  registration  is  much  more  pronounced  in 
this  case,  indicating  that  the  multiscale  hybrid  registration  algorithm  is  particularly 
well-suited  for  registration  problems  involving  large  localized  deformations.  To  vi¬ 
sualize  the  deformation  between  the  images  we  illustrate  the  deformation  computed 
by  the  multiscale  hybrid  registration  algorithm  applied  to  a  grid  image  in  Figure  9. 

Table  5.2  compares  the  correlation  coefficients  between  the  fixed  and  deformed 
images  before  registration,  after  hybrid  multiscale  registration,  after  multiscale  de¬ 
formable  registration,  and  after  landmark  registration.  For  reference,  the  correlation 
coefficient  before  registration  is  0.70. 

5.3.  3D  Registration  of  respiratory  phases.  In  this  section,  we  demonstrate 
that  our  multiscale  hybrid  registration  technique  accurately  registers  three-dimensional 
images.  We  combine  the  80  two-dimensional  images  from  respiratory  phases  1  and 
8  (shown  in  Figure  2)  to  obtain  two  three-dimensional  CT  images,  as  illustrated  in 
Figure  10. 

To  register  the  3D  CT  images  with  one  another,  we  first  extend  the  hierarchical 
multiscale  decomposition  of  [25]  to  3D  images.  Although  the  multiscale  decompo¬ 
sition  presented  in  [25]  was  done  in  two  dimensions  only,  the  hierarchical  multiscale 
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Figure  9.  The  original  and  deformed  grid  image.  The  deformed 
grid  on  the  right  is  obtained  by  applying  the  deformation  computed 
by  the  multiscale  hybrid  registration  algorithm  to  the  grid  image 
on  the  left. 


Registration 

Correlation  Coefficient 

Method 

After  Registration 

Hybrid  multiscale 

0.97 

Multiscale  deformable 

0.93 

Ordinary  deformable 

0.90 

Landmark 

0.84 

Table  1.  The  correlation  coefficients  between  the  fixed  and  de¬ 
formed  images  after  hybrid  multiscale  registration,  after  multiscale 
deformable  registration,  and  after  landmark  registration  for  the 
large  deformation  test  case  shown  in  Figure  7. 


Figure  10.  The  3D  CT  volumes  obtained  upon  combining  the  80 
2D  CT  phase-binned  images  corresponding  to  the  first  (left)  and 
eighth  (right)  respiratory  phases. 


expansion  in  equation  (5)  is  independent  of  the  image  dimensionality.  Upon  decom¬ 
posing  each  3D  image  we  identify  corresponding  landmarks  in  each  image.  For  the 
purposes  of  this  paper  we  identify  the  landmarks  slice-by-slice  (as  in  the  2D  case) 
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and  refer  the  interested  reader  to  [24]  for  an  algorithm  to  automatically  identify 
homologous  regions  in  each  slice.  We  then  register  the  3D  coarse  scales  with  one 
another  using  the  identified  landmarks  and  apply  the  resulting  transformation  to 
the  original  moving  image.  This  coarse-scale  landmark-based  registration  is  imple¬ 
mented  as  a  three-dimensional  registration.  Finally,  we  use  an  iterated  multiscale 
three-dimensional  deformable  B-splines  registration  algorithm  to  complete  the  reg¬ 
istration  process.  The  advantage  of  using  three-dimensional  landmark-based  and 
deformable  registration  procedures  to  register  the  volumes  is  that  three-dimensional 
algorithms  allow  for  transformation  in  three  dimensions  (instead  of  only  allowing 
deformations  restricted  to  a  two-dimensional  plane). 

Figure  11  compares  the  intensity  difference  between  two  sample  slices  before 
and  after  multiscale  hybrid  registration.  A  comparison  of  the  intensity  differences 
before  and  after  registration  demonstrates  that  the  multiscale  hybrid  registration 
method  indeed  recovers  the  difference  between  the  two  images.  Similar  results  were 
obtained  with  all  other  slices.  The  correlation  between  the  images  is  0.8316  before 
registration  and  0.9745  after  registration.  Working  on  a  Dell  Dimension  8400  Intel 
Pentium  4  CPU  (3.40  GHz,  2.00  GB  of  RAM),  the  total  required  computation 
time  for  both  the  3D  multiscale  decomposition  and  multiscale  hybrid  registration 
algorithm  is  on  the  order  of  approximately  20  to  40  minutes,  depending  on  the  data 
set;  this  particular  example  required  29  minutes. 


6.  Robustness  of  the  multiscale  hybrid  registration  algorithm.  In  this  sec¬ 
tion,  we  demonstrate  the  robustness  of  the  multiscale  hybrid  registration  algorithm 
with  several  image  registration  experiments.  All  image  registration  experiments  in 
this  section  are  performed  using  the  fixed  and  deformed  images  shown  in  Figure  7. 
The  fixed  image  in  each  experiment  is  the  image  PIS  16  (first  respiratory  phase,  slice 
16),  and  the  moving  image  is  the  image  obtained  upon  applying  a  large  B-splines 
deformation  to  the  fixed  image.  We  shall  consider  robustness  with  respect  to  the 
location  of  the  landmarks,  number  of  landmarks,  and  several  types  of  noise. 

6.1.  Location  of  landmarks.  In  this  section,  we  demonstrate  the  robustness  of 
the  multiscale  hybrid  registration  technique  with  respect  to  the  location  of  the 
bony  structure  landmarks  identified  using  the  coarse-scale  representations  of  the 
images  to  be  registered.  We  perform  20  image  registration  experiments  in  which 
we  randomly  perturb  the  position  of  the  coarse-scale  landmarks  in  the  fixed  image 
only  by  10  to  20  mm  in  each  trial  (the  image  sizes  are  all  512  x  512).  In  each  trial, 
we  place  the  landmarks  in  the  moving  image  in  exactly  the  same  positions  as  those 
illustrated  in  Figure  4,  but  we  vary  the  location  of  each  corresponding  landmark  in 
the  fixed  image  randomly  by  10  to  20  mm,  as  illustrated  in  Figure  12.  The  moving 
image  landmarks  are  randomly  placed  in  the  blue  circles  shown  in  Figure  12. 

These  experiments  are  designed  to  determine  whether  or  not  the  accuracy  of  the 
hybrid  algorithm  is  dependent  on  precise  matching  between  the  landmarks.  Table 
2  presents  the  mean  square  differences  (MSDs)  and  correlation  coefficients  between 
the  images  after  registration.  The  MSD  between  two  images  A  and  B  is  defined  as 
follows: 


MSD(A,B)  = 

V  2=1 


HYBRID  MULTISCALE  LANDMARK  AND  DEFORMABLE  IMAGE  REGISTRATION  727 


Difference  Difference  After 

Before  Registration  Hybrid  Registration 


Figure  1 1 .  The  difference  between  two  corresponding  slices  of  the 
3D  CT  images  before  and  after  multiscale  hybrid  registration. 


where  N  is  the  total  number  of  pixels  in  each  image,  Ai  is  the  ith  pixel  of  image 
A ,  and  Bi  is  the  ith  pixel  of  image  B.  Note  that  the  optimum  value  of  the  MSD 
is  0,  indicating  exact  matching  between  the  images,  while  poor  matches  between 
the  images  result  in  large  MSD  values.  For  reference,  MSD  before  registration  is 
0.1210  and  the  correlation  coefficient  before  registration  is  0.7924.  For  reference,  we 
also  include  the  MSD  and  correlation  coefficient  after  hybrid  multiscale  registration 
with  exact  correspondence  of  the  landmarks  (i.e.,  no  intentional  perturbation). 

The  results  presented  in  Table  2  demonstrate  that  the  multiscale  hybrid  reg¬ 
istration  technique  is  computationally  robust  with  respect  to  the  location  of  the 
landmarks  used  in  the  coarse-scale  landmark  registration  phase  of  the  algorithm. 
Although  the  precise  location  of  the  landmarks  differs  in  each  of  the  20  trials,  the 
final  registration  results  are  essentially  the  same.  This  is  because  the  coarse-scale 
landmark  registration  is  only  the  first  step  in  the  hybrid  registration  algorithm;  the 
multiscale  deformable  registration  step  fine  tunes  the  result  of  the  landmark  regis¬ 
tration,  thus  correcting  any  anomalies  introduced  in  the  landmark  selection  process. 
This  robustness  indicates  that  the  hybrid  algorithm  is  particularly  well  suited  for 
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Figure  12.  The  perturbed  locations  of  the  coarse-scale  landmarks 
in  the  moving  images.  In  each  trial,  the  landmarks  in  the  moving 
image  are  randomly  placed  in  the  blue  circles. 


clinical  applications.  Although  the  algorithm  does  require  the  identification  of  cor¬ 
responding  landmarks  in  the  coarse-scale  images,  the  algorithm  is  not  sensitive  to 
the  exact  location  of  the  landmarks. 

This  robustness  to  the  location  of  the  landmarks  is  one  of  the  most  notable 
features  of  our  hybrid  multiscale  landmark  and  deformable  registration  algorithm. 
Most  landmark-based  registration  algorithms  require  precise  identification  of  cor¬ 
responding  landmarks  and  are  thus  tedious,  time-consuming,  and  difficult  to  im¬ 
plement,  even  for  a  medical  expert.  However,  since  the  accuracy  of  our  hybrid 
algorithm  is  not  dependent  on  precise  correspondence  of  the  landmarks  in  the  fixed 
and  moving  images,  the  identification  of  the  landmarks  can  be  completed  quickly 
or  implemented  automatically. 

6.2.  Number  of  landmarks.  In  this  section,  we  consider  the  effect  of  the  number 
of  pairs  of  landmarks  used  in  the  landmark-based  registration  step  on  the  accuracy 
and  efficiency  of  the  multiscale  hybrid  registration  algorithm.  Table  3  presents  the 
correlation  coefficients  between  the  images  after  registration  and  total  registration 
time  (including  the  multiscale  decomposition,  landmark-based  registration  of  the 
coarse  scales,  and  deformable  registration)  for  image  registration  experiments  using 
a  varying  number  of  landmarks.  The  locations  of  the  landmarks  are  illustrated  in 
Figure  13.  The  results  presented  in  Table  3  indicate  that  accurate  registration  can 
be  obtained  using  four  pairs  of  corresponding  landmarks. 

6.3.  Noise.  In  this  section,  we  demonstrate  the  robustness  of  the  multiscale  hybrid 
registration  algorithm  with  respect  to  the  presence  of  noise  in  the  images  to  be 
registered.  Our  multiscale  registration  algorithms  [15],  [16]  were  initially  developed 
in  the  context  of  registration  of  noisy  medical  images.  We  consider  three  types  of 
noise  models:  Gaussian,  multiplicative,  and  impulse  noise. 

6.3.1.  Gaussian  noise.  Gaussian  noise  is  an  independent  additive  noise  model  in 
which  the  observed  (noisy)  image  f{x)  is  the  sum  of  the  true  image  s(x)  and  the 
noise  n(x): 
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Trial  Number 

MSD  After 
Registration 

Correlation  Coefficient 
After  Registration 

Exact  Landmark 

0.0056 

0.9686 

Correspondence 

1 

0.0055 

0.9643 

2 

0.0054 

0.9892 

3 

0.0059 

0.9741 

4 

0.0048 

0.9823 

5 

0.0050 

0.9696 

6 

0.0047 

0.9728 

7 

0.0049 

0.9781 

8 

0.0051 

0.9690 

9 

0.0055 

0.9704 

10 

0.0053 

0.9668 

11 

0.0055 

0.9712 

12 

0.0049 

0.9842 

13 

0.0047 

0.9901 

14 

0.0050 

0.9834 

15 

0.0052 

0.9745 

16 

0.0051 

0.9736 

17 

0.0054 

0.9963 

18 

0.0047 

0.9896 

19 

0.0049 

0.9800 

20 

0.0048 

0.9773 

Table  2.  The  MSDs  and  correlation  coefficients  between  the  fixed 
and  deformed  images  after  multiscale  hybrid  registration.  Each 
trial  number  represents  a  different  perturbation  of  the  landmarks 
in  the  fixed  image. 


Figure  13.  The  locations  of  the  landmarks  used  in  the  experi¬ 
ments  presented  in  Table  3. 
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Number  of  Pairs 
of  Landmarks 

Correlation  Coefficient 
After  Registration 

Total 

Computation  Time 

2 

0.9365 

99  seconds 

3 

0.9488 

89  seconds 

4 

0.9686 

75  seconds 

5 

0.9598 

77  seconds 

6 

0.9486 

85  seconds 

7 

0.9565 

86  seconds 

Table  3.  The  correlation  coefficients  between  the  fixed  and  de¬ 
formed  images  after  hybrid  registration  and  total  computation  time 
required  for  hybrid  registration  for  image  registration  experiments 
using  a  varying  number  of  landmarks. 


Figure  14.  The  fixed  and  deformed  images  with  Gaussian  noise 
of  mean  0  and  variance  0.1. 


f(x)  =  s(x) +n(x),  (6) 

where  n(x)  is  uniformly  distributed  random  noise  with  mean  (i  and  variance  v. 
Figure  14  illustrates  the  fixed  and  moving  images  with  added  Gaussian  noise  of 
mean  0  and  variance  0.1 

Table  4  presents  the  correlation  coefficients  between  the  noisy  images  shown 
in  Figure  14  before  registration,  after  ordinary  deformable  registration,  and  after 
hybrid  registration. 

6.3.2.  Multiplicative  noise.  Multiplicative  noise,  or  speckle  noise,  is  commonly  ob¬ 
served  in  medical  imaging.  It  is  defined  by  the  following  model,  where  s(x)  denotes 
the  actual  image  and  f(x)  denotes  the  observed  (noisy)  image: 

f{x)  =  s(x)  +r)(0,5)  ■  s(x),  (7) 

where  rj(0,6)  is  uniformly  distributed  random  noise  of  mean  0  and  variance  5. 
where  n(x)  is  uniformly  distributed  random  noise  with  mean  (i  and  variance  v. 
Figure  15  illustrates  the  fixed  and  moving  images  with  added  multiplicative  noise 
of  mean  0  and  variance  0.2 
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Gaussian  Noise 

Correlation  Coefficient 

Computation  Time 

Before  registration 

0.1830 

Deformable  registration 

0.2057 

304  seconds 

Multiscale  deformable  registration 

0.8915 

92  seconds 

After  hybrid  registration 

0.9149 

81  seconds 

Table  4.  The  correlation  coefficients  between  the  noisy  fixed  and 
moving  images  shown  in  Figure  14  before  registration,  after  ordi¬ 
nary  deformable  registration,  and  after  multiscale  hybrid  registra¬ 
tion. 


Figure  15.  The  fixed  and  deformed  images  with  multiplicative 
noise  of  mean  0  and  variance  0.2. 


Multiplicative  Noise 

Correlation  Coefficient 

Computation  Time 

Before  registration 

0.3527 

Deformable  registration 

0.4683 

295  seconds 

Multiscale  deformable  registration 

0.9011 

89  seconds 

Hybrid  multiscale  registration 

0.9366 

76  seconds 

Table  5.  The  correlation  coefficients  between  the  noisy  fixed  and 
moving  images  shown  in  Figure  15  before  registration,  after  ordi¬ 
nary  deformable  registration,  and  after  multiscale  hybrid  registra¬ 
tion. 


Table  5  we  presents  the  correlation  coefficients  between  the  noisy  images  shown 
in  Figure  15  before  registration,  after  ordinary  deformable  registration,  and  after 
hybrid  registration. 

6.3.3.  Impulse  noise.  Impulse  noise,  or  salt  and  pepper  noise,  is  noise  that  resembles 
salt  and  pepper  granules  randomly  distributed  over  the  image.  Impulse  noise  is 
typically  defined  by  the  following  model,  where  s(x)  denotes  the  actual  image  and 
f(x)  denotes  the  observed  (noisy)  image: 
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Figure  16.  The  fixed  and  deformed  images  with  impulse  noise  of 
density  0.2. 


Impulse  Noise 

Correlation  Coefficient 

Computation  Time 

Before  registration 

0.2226 

Deformable  registration 

0.2451 

333  seconds 

Multiscale  deformable  registration 

0.9042 

105  seconds 

Hybrid  registration 

0.9377 

82  seconds 

Table  6.  The  correlation  coefficients  between  the  noisy  fixed  and 
moving  images  shown  in  Figure  16  before  registration,  after  ordi¬ 
nary  deformable  registration,  and  after  multiscale  hybrid  registra¬ 
tion. 


f(x) 


s(x),  with  probability  1  —  5, 

r)(x),  with  probability  5, 


(8) 


where  rj(x)  is  an  identically  distributed,  independent  random  process  which  sets 
corrupted  pixels  alternatively  to  zero  (black)  or  one  (white);  unaffected  pixels  re¬ 
main  unchanged.  An  arbitrary  pixel  is  affected  by  noise  with  probability  5,  and  not 
affected  with  probability  1  —  5.  We  refer  to  S  as  the  impulse  noise  density,  since 
adding  impulse  noise  of  density  S  to  an  image  f{pc)  affects  approximately  S  •  siz e(/) 
pixels.  Figure  16  illustrates  the  fixed  and  moving  images  with  added  impulse  noise 
of  density  0.2 

Table  6  presents  the  correlation  coefficients  between  the  noisy  images  shown  in 
Figure  16  before  registration,  after  ordinary  deformable  registration,  after  multi¬ 
scale  deformable  registration,  and  after  hybrid  registration,  as  well  as  the  total 
computation  time  required  for  each  registration. 


6.4.  Discussion.  The  results  presented  in  this  section  demonstrate  the  accuracy, 
efficiency,  and  robustness  of  the  hybrid  multiscale  landmark  and  deformable  reg¬ 
istration  technique  for  registration  of  medical  images  in  both  two  and  three  di¬ 
mensions.  In  particular,  the  hybrid  technique  is  shown  to  be  more  accurate  than 
ordinary  landmark  registration,  ordinary  deformable  registration,  and  multiscale 
deformable  registration,  especially  in  the  case  of  registering  images  that  contain 
large  localized  deformations  and/or  noise.  Because  the  accuracy  of  the  hybrid 
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technique  is  not  sensitive  to  precise  identification  of  the  location  of  correspond¬ 
ing  landmarks,  the  hybrid  registration  method  is  a  significant  improvement  over 
ordinary  landmark-based  registration  methods,  which  are  known  to  be  dependent 
on  exact  spatial  correspondence  of  landmarks.  Additionally,  the  hybrid  technique 
accurately  registers  deformed  images  that  contain  noise  levels  large  enough  that 
ordinary  registration  techniques  fail  to  produce  meaningful  results. 

The  development  of  our  hybrid  multiscale  technique  was  motivated  by  the  ob¬ 
servation  that  the  correspondence  between  some  of  the  regions  in  the  moving  image 
(such  as  bony  structures)  can  be  easily  identified  visually,  while  the  correspondence 
between  other  regions  (such  as  those  that  contain  tissue  deformation,  breathing 
movement,  or  lack  of  distinct  image  features)  is  less  obvious.  Thus,  rather  than 
approach  the  mapping  of  the  two  types  of  regions  equally,  as  with  ordinary  image 
registration  algorithms,  mapping  should  be  approached  separately  for  each  region. 
Our  proposed  technique  combines  landmark  registration  and  deformable  registra¬ 
tion  so  that  prior  knowledge  about  the  correspondence  between  the  images,  such 
as  visual  identification  of  corresponding  landmarks,  can  be  incorporated  into  the 
deformable  registration  process.  The  use  of  the  coarse  scales  (which  contain  only 
the  main  shapes  and  general  features  of  the  images)  obtained  via  the  hierarchical 
multiscale  image  decomposition  of  [25]  enables  quick  and  easy  identification  of  cor¬ 
responding  landmarks,  and  the  deformable  registration  component  of  the  algorithm 
fine-tunes  the  registration  result  produced  by  the  coarse-scale  landmark  registration. 

7.  Conclusions.  In  this  paper,  we  proposed  a  novel  hybrid  two-step  image  reg¬ 
istration  technique  combining  the  hierarchical  multiscale  image  decomposition  of 
[25]  with  landmark  and  deformable  registration  methods.  The  hybrid  approach  has 
several  major  advantages  over  ordinary  registration  algorithms.  It  allows  the  practi¬ 
tioner  to  incorporate  a  priori  knowledge  of  corresponding  bony  or  other  anatomical 
structures  into  the  registration  process  using  a  coarse-scale  landmark  registration. 
Because  the  coarse  scales  of  an  image  contain  only  its  main  shapes,  the  coarse-scale 
landmark  registration  phase  of  the  hybrid  algorithm  is  easy  to  implement  and  is 
computationally  efficient.  The  transformation  produced  by  the  coarse-scale  land¬ 
mark  registration  is  used  as  the  starting  point  for  a  deformable  registration,  which 
significantly  decreases  the  computation  time  required  for  convergence  of  the  de¬ 
formable  registration  algorithm.  The  hybrid  method  was  applied  to  both  two-  and 
three-dimensional  deformable  image  registration  problems,  and  our  results  demon¬ 
strated  that  the  hybrid  registration  is  more  accurate  that  ordinary  landmark  and 
deformable  registration  methods,  and  that  the  technique  significantly  improves  the 
accuracy  of  registration  of  images  that  contain  large  localized  deformations.  The 
hybrid  algorithm  is  very  robust  with  respect  to  the  location  of  the  landmarks  used 
in  the  coarse-scale  landmark  registration  phase  of  the  algorithm,  indicating  that 
the  accuracy  of  the  technique  is  not  dependent  on  precise  spatial  correspondence  of 
landmarks.  Compared  to  ordinary  landmark  registration,  this  method  significantly 
reduces  the  time  required  for  the  difficult  and  tedious  task  of  landmark  selection. 
The  hybrid  method  is  also  robust  with  respect  to  the  presence  of  several  different 
types  of  noise,  and  greatly  improves  the  accuracy  of  registration  of  images  that 
contain  noise  or  other  artifacts. 

Appendix  A.  Implementation  of  the  multiscale  decomposition.  As  de¬ 
scribed  in  [25],  the  initial  scale  Ao  should  capture  the  smallest  oscillatory  scale 
in  /,  given  by 


734 


D.  PAQUIN,  D.  LEVY,  AND  L.  XING 


|w-i>o°  —  ” — ’ 


< 


2A0  "  -  A0’ 

where  W~1,oc  is  the  Sobolev  space  with  norm  given  by: 


(9) 


W-1’°° 


sup 

9 


/ 


f(x)g(x) 


dx 


|  1^1  |  W1’1 

where  ||g||vpi,i  tm  ||V<7||li.  However,  in  practice  we  may  not  be  able  to  determine 
the  size  of  \\f\\w-i,oo,  so  we  determine  the  initial  choice  of  Ao  experimentally.  Fol¬ 
lowing  [25],  for  the  applications  presented  in  this  paper,  we  will  use  Ao  =  0.01  and 
\j  =  A02J. 

We  follow  the  numerical  algorithm  of  [25]  for  the  construction  of  our  hierarchical 
decomposition.  In  each  step  we  use  finite- difference  discretization  of  the  Euler- 
Lagrange  equations  associated  with  the  J(vj,Xj^i)  to  obtain  the  next  term,  Uj+ 1, 
in  the  decomposition  of  the  image  /.  Due  to  the  singularity  when  \Vu\\  =  0,  we 
replace  J(/,  A)  by  the  regularized  functional 


Je(/,A):=  inf 

U+V  =  f 


bill2  + 


A 


e2  +  |  V^|2  dx  dy 


}■ 


(10) 


and  at  each  step,  we  find  the  minimizer  u\  of  Je.  The  Euler-Lagrange  equation  for 

■I'd,  A)  is 


n,  -  -div 


Vux 


ye^iv^ib 

with  the  Neumann  boundary  conditions: 


=  /  in  ft  , 


dux 

dn 


=  0, 


(11) 


an 


where  dfl  is  the  boundary  of  the  domain  Q  and  n  is  the  unit  outward  normal.  We 

k 

thus  obtain  an  expansion  /  ~  uji  where  the  Uj  are  constructed  as  approximate 

3=0 

solutions  of  the  recursive  relation  given  by  the  following  elliptic  PDE: 


uj+ 1 


-div 


\7u 


'3  + 1 


VlG 


.  _ _  = — =-div  -==  (12) 

2^+i  v^  +  iWi+iiy  2^i  We2  +  iv«iiv 

To  numerically  implement  the  method,  we  cover  the  domain  Q  with  a  grid  (xi  :  = 
ih,yj  :=  jh),  and  discretize  the  elliptic  PDE  of  Eq.  (12)  as  follows:. 


^i,j  fi,j  T 

1 

2 h? 


"i,3 


di i  j  Ui— i 


d/e2  +  ( D^~u~/p  +  (Doj/Wi.j)2  \/e2  +  ( U~u~ )2  +  (£>0j/Wi- i,j)2 

2^2  \A2  +  +  (D+^tijj)2  \/e2  +  (-fV^ij-i)2  +  (D-yWjj)2 


(13) 
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where  D+,  D_,  and  Dq  denote  the  forward,  backward,  and  centered  divided  dif¬ 
ferences,  respectively.  To  solve  the  discrete  regularized  Euler-Lagrange  equations 
(13),  we  use  the  Gauss-Siedel  iterative  method  to  obtain: 


.n+l 


fij  + 


1 

2 h? 

1 

2 h? 


ui+l,j  Ui,j 


UlV  Ui- 1,3 


_\f<?  +  (D+xU?tj)2  +  (DoyW^)2  /e 


+  {D_xuljy  +  {DQyuU,j)\ 


_  7/n+l 

Ui,j+ 1  Ui,j 


7/n+!  _  o,n 

ui,j  ai,j- 1 


V*2  +  (A>*«&)2  +  (-O+s/M^)2  y/*2  +  (Ate^-i)2  +  T-^)2. 


(14) 


To  satisfy  the  Neumann  boundary  conditions  (11),  we  first  reflect  /  outside  11  by 
adding  grid  lines  on  all  sides  of  Q.  As  the  initial  condition,  we  set  u\  -  =  fij.  We 
iterate  this  numerical  scheme  for  n  =  0,1 , ...  TV  until  ||Rno°  —  Rn°°_1||  is  less  than 
some  preassigned  value  so  that  is  an  accurate  approximation  of  the  fixed  point 
steady  solution  u\. 

Finally,  we  denote  the  final  solution  u\  :=  To  obtain  the  hierarchical 

multiscale  decomposition,  we  reiterate  this  process,  each  time  updating  /  and  A  in 
the  following  way: 


fn 
A  n 


/ current 
2  A  current  • 


(15) 


That  is,  at  each  step,  we  apply  the  J (/current  2 A)  minimization  to  the  residual 
/current  ~  ux  of  the  previous  step.  Taking  A j  =  Ao2J,  we  obtain  after  k  steps  a 
hierarchical  multiscale  decomposition  /  =  u\0  +  u\1  +  . . .  +  u\k  +  v\k ,  where  we 
write  u\j  =  Uj.  We  call  the  Uj,  j  =  1,  2, . . . ,  k  the  components  o/  /,  and  the  vj~  the 
residuals.  For  ease  of  notation,  given  an  image  /,  we  let  Ck(/)  denote  the  kth  scale 
o/  the  image  /,  k  =  1, . . . ,  m: 


k- 1 

W)  =  5>(/).  (16) 

i= 0 

Thus  Ck(A)  will  denote  the  kth  scale  of  the  image  A ,  and  Ck(B)  will  denote  the 
kth  scale  of  image  B. 


A.l.  Three-dimensional  implementation.  To  implement  the  iterated  multi¬ 
scale  decomposition  in  three  dimensions,  we  cover  the  image  domain  Q  with  a 
grid  (xi  :=  ih,yj  :=  jh,Zk  :  =  kh ),  and  let  D+,  D_,  and  Dq  denote  the  forward, 
backward,  and  centered  divided  differences,  respectively.  Then  the  3D  extension  of 
the  iterated  multiscale  decomposition  given  by  equation  (14)  in  Section  A  is: 
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